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We study the reconstruction of visual stimuli from spike trains, 
recording simultaneously from the two HI neurons located in the 
lobula plate of the fly Chrysomya megacephala. The fly views two 
types of stimuli, corresponding to rotational and translational dis- 
placements. If the reconstructed stimulus is to be represented by 
a Volterra series and correlations between spikes are to be taken 
into account, flrst order expansions are insuflicient and we have 
to go to second order, at least. In this case higher order corre- 
lation functions have to be manipulated, whose size may become 
prohibitively large. We therefore develop a Gaussian-like repre- 
sentation for fourth order correlation functions, which works ex- 
ceedingly well in the case of the fly. The reconstructions using this 
Gaussian-like representation are very similar to the reconstructions 
using the experimental correlation functions. The overall contri- 
bution to rotational stimulus reconstruction of the second order 
kernels - measured by a chi-squared averaged over the whole ex- 
periment - is only about 8% of the flrst order contribution. Yet if 
we introduce an instant-dependent chi-square to measure the con- 
tribution of second order kernels at special events, we observe an 
up to 100% improvement. As may be expected, for translational 
stimuli the reconstructions are rather poor. The Gaussian-like rep- 
resentation could be a valuable aid in population coding with large 
number of neurons. 



1 Introduction 



Living animals have to reconstruct a representation of the external world 
from the output of their sensory systems in order to correctly react to the 
demands of a rapidly varying environment. In many cases this sensory output 
is encoded into a sequence of identical action potentials, called spikes. If we 
represent the external world by a time-dependent stimulus function s(t), the 
animal has to reconstruct s{t) from a set of spikes. This decoding procedure 
generates an estimate Se{t) of the stimulus like a digital-to-analog converter. 



Figure 1: Motion sensitivity of the two HI neurons. Eacli eye sees a monitor 
displaying a rigidly moving bar pattern. The stimuli in this figure correspond to 
a translational motion in which both neurons are excited. Inverting the stimulus 
shown by monitor Ml would generate a rotational stimulus, which now inhibits 
the response of the left neuron. Electrodes record extracellularly from each HI. 

Here we study this decoding procedure in a prominent example of spiking 
neurons: the two HI neurons of the fly Chrysomya megacephala. The fly has 
two compound eyes with their associated neural processing systems (Hansen, 
1981, 1982, 1984). Motion detection starts at the photoreceptor cells, eight 
of them located in each one of the ~ 5000 ommatidia of each compound 
eye. They effect the transduction of photons into electrical signals, which are 
propagated via the lamina and medulla to the lobula plate. This neuropil is - 
inter alia - composed of horizontally and vertically directionally sensitive wide 
field neurons. The HI neurons are horizontally sensitive and are excited by 
ipsilateral back to front motion and inhibited by oppositely moving stimuli. 
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Rotational and Translational Rasters 
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Figure 2: Raster plot for the two HI neurons, showing their complementary ac- 
tion under rotational and translational stimuli. The same time-dependent stim- 
ulus s(t) is repeatedly shown to the fly, the horizontal time axis running from 
time zero to 5000 bins = 10 seconds and the vertical axis showing the repetition 
number. The responses of the neurons are shown as a raster, where each dot 
represents a spike. The right HI sees a stimulus Sr{t) and the left one sees si(t). 
Rotational stimuli s,.(t) = si{t) = s{t): (Rl) spikes from right HI and (R2) 
spikes from left HI. Translational stimuli Sr{t) = —si{t) = s{t): (Tl) spikes 
from right HI and (T2) spikes from left HI. Inset to (R2): Spikes from right 
HI, fly subjected to sign reversed stimuli in order to simulate raster (R2). 

Each HI neuron projects its axon to the contralateral lobula plate, exciting 
there two horizontal and two centrifugal cells. These cells mediate mutual 
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inhibition between the two HI neurons (Haag, Vermeulen, & Borst, 1999; 
Haag & Borst, 2001; Farrow, Haag, & Borst, 2003; Haag & Borst, 2008; 
Krapp, 2009)^. We subject the fly to rotational and translational stimuh - see 
Figure 1. If the fly rotates around a vertical axis, say clockwise when looking 
down the axis, the left neuron is inhibited and the right one is exited, so 
that the two neurons become an efficient rotational detector (Hansen, 1984). 
This can be seen in Figure 2 (Rl) & (R2). Even when recording only from 
the ipsilateral HI, one can simulate the response of the contralateral HI. In 
fact, since the two HI cells have mirror symmetric directional sensitivities, 
the sign flipped stimulus induces a response in the ipsilateral HI typical for 
the contralateral HI cell (Rieke, Warland, Steveninck, & Bialek, 1997). The 
inset in (R2) shows this to be true to a very good approximation. 

In forward translation none of HI neurons is excited, corresponding to 
the low spike density regions in the raster-plots of Figure 2 (Tl) & (T2). 
In backward translation, both Hi's are excited and we expect a strong inhi- 
bition. Yet the spike rate is comparable to rotational excitation - compare 
Figure 2 (Rl) & (Tl). Numerical computation conflrms this visual impres- 
sion. Nevertheless in translation the two Hi's flre mainly in sync, which 
leads to subtle differences with respect to rotation. As a consequence, our 
reconstructions will be much poorer for the translational case - see section 5. 

If we want to take correlations between spikes into account, instead of 
treating them independently, we have to go at least to second order stimulus 
reconstructions. These require the computation of higher order spike-spike 
correlation functions and a subsequent matrix inversion. If one records from 
many neurons simultaneously, the size of these matrices may soon become 
prohibitively large. Here we present an efficient representation of these higher 
order correlation functions in terms of second order ones. The reconstruc- 
tion now costs far less computationally, avoids large matrix inversions and 
gives excellent results. We test the quality of our reconstructions under both 
rotational and translational stimuli. 

If this representation holds more generally, it may well make population 
coding computationally more tractable. We briefly discuss a perturbation 
scheme, which allows a stepwise inclusion of small effects. 

^Although experimental work has focusscd on the vertical system, one expects analog 
results for the horizontal one. 
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2 Stimulus reconstruction from spike trains 



Suppose we want to reconstruct the stimulus from the response of a single 
HI neuron. We represent this response as a spike train p{t) = ^^J^i S{t — ti), 
which is a sum of delta functions at the spike times ti. Ng is the total number 
of spikes generated by the neuron during the experiment. 

The simplest reconstruction extracts the stimulus estimate via a linear 
transformation, see e.g. (Rieke et al., 1997; Bialek, Rieke, Steveninck, & 
Warland, 1991), 

/oo 
hiT)pit-T)dT, (2.1) 
-oo 

with the kernel ki{t) to be determined. 

For simplicity we effect an acausal reconstruction, i.e. we integrate from 
— oo to +CXO. Essentially the same results are obtained in a causal recon- 
struction. One way to implement causality proceeds to estimate the stim- 
ulus at time t, using as input the spike train up to time t + to. For the 
fly to has to be > to 30 milliseconds. In this case equation 2.1 would read: 

Se{t) = f^^^ki{T)p{t-T)dT. 

Equation 2.1 is the first term of a Volterra series (Martin, 2006): 

roc poo 

Se{t)= ki{T)p{t-r)dT+ k2{ri,r2)p{t-n)p{t-T2)dridT2 + ... {2.2) 



There is no convergence proof for this expansion, but heuristically we may 
say that it should be a valid approximation, if the average number of spikes 
per correlation time Tc, 

V = {r)rc, (2.3) 

is small (Rieke et al., 1997). Here (r) is the mean spike rate and Tc a typical 
signal correlation time. For small t] each spike gives independent information 
about the stimulus. In our case t] ~ 0.6 — 0.8, which is of the order of unity, 
so that higher order effects might be relevant. 

The first order term, being proportional to Ylf" ~ ^i)' independently 
adds contributions for each spike. Yet it is well established that pairs of spikes 
carry a significant amount of additional information beyond the single spike 
contributions (Brenner, Strong, Koberle, Bialek, & Steveninck, 2000). This 
motivates the addition of the second order kernel k2{Ti,T2), which includes 
correlations between up to two spikes. 
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In order to obtain the kernels ki and k2 we choose to minimize the fol- 
lowing functional - the x*^^"* error - 

^'\k„k,) = {[dt[s,{t)-s{t)]'). (2.4) 



The brackets stand for an ensemble average with respect to the distribution 
of all possible stimuli in a given experiment. In a long experiment we average 
over ~ 10^ time windows of size T^. Typically Ty^ ~ 100 milliseconds 
- see section 7 for details. For ease of presentation, in the following our 
discussions will always refer to the rotational setup, unless explicitly stated 
otherwise as in section 5. 

Since the functional 2.4 is quadratic, the equations minimizing x^'^H^i, ^2) 

dx'^'ydk,=0,j = l,2 (2.5) 

are linear in the unknowns ki, k2- E.g., if we keep only ki, using therefore 
equation 2.1, we get: 

~ (g(^)*p(^)) ,^ „^ 

{p{uj)*p{uj)) 

where Fourier transforms are defined as F{u!) = J dtF{t)e^^^. 

We may include the second order term k2-, either as a correction to the 
first order reconstruction Si{t) = ki-k pit) ^, or one may solve the coupled 
system 2.5. 

If we record simultaneously from left and right HI, we obtain two spike 
trains pi{t) and p2{'t). The expansion equation 2.2 generalizes to 

Se{t) = Ki^pi{t)+K2^P2{t) + 

Kn -k pi -k pi{t) + Ku-k pi-k p2{t) + K22 -k P2 k p2it) + .... (2.7) 

Here we have included the kernel K12, which encodes effects correlating pi 
and P2 ^. Notice that K12 = K21. 

To first order, keeping only Ki and K2 in the expansion 2.7, we get the 
following equations: 

2 

SR^iuj) = J2Kb{u^)Rab{u^), a = 1,2 (2.8) 
6=1 



^The symbol * stands for a convolution as in equation 2.1. 

^Notice that we have not orthogonalized our expansion equation 2.2. so that there are 
^11(^1,^1) terms, which could have been absorbed in Ki{t) and similarly for K2{t). 
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where 

ma{uj) = I dtdt'{s{t')pa{t' - t))e'^\ (2.9) 

and 

Rab{tiM) = I dt{pa{t - ti)pbit - ta)), a, 6 = 1, 2. (2.10) 



Due to time-translation invariance Rab{ti,t2) is only a function of the dif- 
ference: Rabiti,t2) = Rabiti — ^2) and Rabijjj) = J dtRab{t)e^^ . Analogous 
properties hold for all the following correlation functions involving only p{t). 
The solution of equations 2.8 yields 

ka{uj) = {La{uj)Raa - La{u) Raa{cu)) / A, a = 1,2 (2.11) 

where 

Laiuj) = {siuj)p:{iu)), A = RnR22 - R12R21 (2.12) 
and d = 3 — a. We obtain the first order reconstruction as 

s,{t) = K^^p,{t) + K2^P2{t). (2.13) 

Since the second order contribution turns out to be small, we treat it 
as a perturbation to the first order reconstruction. We therefore expand 
S2{t) = s{t) — Si(t) as: 

S2{t) = Kii-k pi-k pi{t) + Ki2-k pl-k p2{t) + K22'k P2'k p2{t). (2.14) 

We now have to solve the following equations 

3i'ab{tl,t2) = / dtsdU Kcd{tl,t2)R'^la{ti,t2,t3,U), (2.15) 



c4=l 



where 



3i^bitl,t2)= j dt{s2{t)pa{t-t^)ph{t-t2)), (2.16) 
R''lditl,t2,h,U) = j dt{pa{t-h)pi,{t-t2)pc{t-h)pd{t-h)). (2.17) 

Although the system 2.15 is linear, the matrices to be inverted may be 
very large. We have to invert the matrix 

M^jS^ ^ R^ScditiMMM). (2.18) 
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where A,B are compound indices A = [ab], B = [cd] labeling the neurons. 
T = [^1,^2],^' = [^35^4] are compound time indices of size each. If we 
compute the correlation functions using a time window of Ty^ = 128 bins, with 
binsize = 2 milliseconds, then the size of M^§' is ~ 128^ x 2"^ ~ 5 x 10^. 
The matrices to be inverted may become prohibitively large, especially if we 
record from more than just two neurons ^. 

We therefore present below a Gaussian-like representation of -R^^j,^ with 
a small number of parameters and which requires no large matrix inversion. 

3 Gaussian-like (Gl) representation for 4-point 
functions 

In this section we present a representation of the 4-point function i?^^cd 

(2) 

terms of the 2-point function -R),^ , which is surprisingly good and which 
avoids the computation of the large matrices 2.18. 

If our spike-generating process were Gaussian, we would have the follow- 
ing structure for R^'^^: 

R^^\l, 2, 3, 4) = R{1, 2)R{3, 4) + R{1, 3)R{2, 4) + R{1, A)R{2, 3) 

-2(p(t))^ (3.19) 

where {p(t)) is just a constant, due to time-translation invariance^. 
This suggests the following representation for R^^^: 

(1, 2, 3,4) = A [R{1, 2)i?(3, 4) + R{1, 3)R{2, 4) + 

i?(l,4)i?(2,3)] - B, (3.20) 

where A and B are constants to be adjusted^. 
For two neurons we get the representation: 

Rabcdil, 2, 3, 4) = [Rab{l, 2)R,d{3, 4) + i?„,(l, 3)Rm{2, 4) + 

Radi^, 4)i?bc(2, 3)]Aabcd + Babcd (3-21) 

^We may solve the above system in Fomier space and select a subset of frequencies in 
order to reduce the size of the system. 

^We write (1, 2, . . .) instead of (ti, t2, ■ • ■)• 

^Any structure built only from could be used for our method to work. 
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Figure 3: Window-size dependence of parameters Ann and -Bun. Similar 
behavior is found for the other parameters A and B. Notice that variations are 
on the 0.05% level. 



with a, b,c,d = 1,2 and Aabcd, Babcd constants to be determined. 

The usefulness of our Gl-representation scheme depends on the quahty 
of the 4-point functions obtained, which in turn hinges on the knowledge of 
the constants Aabcd and Babcd- There would be no point, if this required the 
computation of 4-point functions in large window sizes and a fitting procedure 
using these windows - exactly what we wanted to avoid. We therefore fit the 
constants Aabcd and Babcd for a sequence of window sizes T^, ranging from 
10 to 128 bins, using -Riiii(ti,t2 = ^3 = ^4 = 1) to fit to the experimental 
data. As can be seen in Figure 3, at least in the fly's case, the dependence 
of the parameters Aabcd) Babcd on is only 0.05% and therefore completely 
negligible. The constants Aabcd and Babcd can therefore be computed very fast 
in small windows. In Figure 4 we plot the fits to the first row -Riiii(ti,t2 = 
^3 = ^4 = 1) and its Gl approximation. As advertised we obtain a perfect fit. 
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Figure 4: 4-point functions and its Gl approximation. We plot -Riiii(ti,t2 = 
^3 = ^4 = 1) for window-size = 64 bins. The black continuous line is 
the experimental 4-point function. The dashed line is its Gaussian approxima- 
tion without parametrization using equation 3.19. The circles represent its Gl 
approximation 3.20. 



In Figure 5 we show the Gl approximation for the Ruu{ti, t2, ^3 = ^4=1 
and its experimental version, which emphasizes the quality of the approx- 
imation. Using the same parameters for the other entries of -Rim and for 
i?2222 results in a fitting error about 20 % larger. 

One of the utilities of this representation will become apparent, once we 
deal with the solution of equation 2.15 in the next section. 

4 A convenient set of functions to solve for 
second order kernels 

At this point it is convenient to introduce a complete set of basis functions 
= l,2,..,nf to expand our variables in. We thus trade continuous 
time-arguments for discreet Greek indices. We expand our second order 
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Figure 5: (A) Gl approximation and (B) experimental 4-point function for 

^2, ^3 = = 1) for window-size = 64 bins. 



kernels as: 



,ab 



fl.U 

We also expand our correlation functions: 

3i^at!ih,h) =J2^t fMMt2) 



(4.22) 



(4.23) 



and 



-^ifeL(^l'^2,^3,^4) 
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Figure 6: 2-point correlation functions rii(t) = (pi(ti — t)pi(ti)), r22(t) = 
(P2(^i — t)p2iti)) and ri2(t) = (pi(ti — t)p2(^i))- The central peak is absent in 
the mixed correlation function ri2(t). 



In order to efficiently compute our second order kernels it is crucial to 
select an adequate set for f^(t), fi = 1, 2, .., n/. 

Depending on the case, it may be sufficient to use a small number n/ 
of functions /^(t) to get a useful representation. If nj has only a slight 
dependence on window size T^, this would allow one to increase without 
further computational costs. 

Often a Fourier expansion is used, i.e. = e*"^*. But we may exploit 
our liberty to choose the functions in a more profitable way. Since our 2- 
point function R(ti,t2) is real, positive^ and symmetric in ti,t2, it posses a 

''in case this is not true, we just add a convenient constant. 
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complete set of eigenf unctions hfj^it): 

j dt2R{ti,t2)h^{t2) = r^h^iU) (4.25) 

with eigenvalues r^, = 1, . . . , A^^. We now choose our functions as f^if) = 
hfj_(t) / y/rj^, which satisfy: 

j dhdt2Uti)R{ti,t2)Mt2) = 5^,. (4.26) 

This choice will avoid large matrix inversions, if at least part of our higher 
order correlation functions can be built from R{ti,t2). 

Substituting the expansions 4.23 and 4.24 into equations 2.15, we get a 
linear system to be solved for PJJ^: 

cd,al3 

In order to avoid cluttering our expressions with indices, we introduce 
our representation first for one neuron only, suppressing thus the indices 
a,b, .., all set to 1. We choose our functions ffi{t) to diagonalize R^^{ti,t2) = 

{Piih)p2it2)y. 

J dhdt2Uh)R'\h,t2)Ut2) = 5^,. (4.28) 

The first of equations 3.21 for TZj^J^/^ becomes 

'n^,^ap = M^^i•^^a^3 + 25^a5^/3) - 2B Uanpn^n^, (4.29) 

where = / dtf^,{t){p{t)). 

Using this expression and the shorthand S^^ = Sj^l in equations 4.27, we 
get the following equations for the unknown coefficients V^^^, = Vj^]^ 

S^, = A[tr{V)5^, + 2V^,] - 2BDnn n,n^, (4.30) 

where tr{V) = Y.^V^,^ and D„„ = Y.^pnaVo,pni3. The sums over /i, 
run from 1 to bins. 

This system can now easily be solved by: 

1. taking the trace over pv to compute tr{T>) = D and 
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2. multiplying by n^^Uy to compute -D„„. 
We get 



with 



where 



V^,, = [S,,JA -D5^, + 2Bn^n, D J/2, (4.31) 

D = [2(1 - ^4)5^;. + 2n2 n^S^,n,]/A, (4.32) 
Dnn = [{n + 2)n^S^unu - S^^n2]/A, (4.33) 



A = 2(T„ + 2)(l-n4)+2n2,n2 = J]n^n^,n4 = (ns)'. (4.34) 

For two neurons we now have to decorate our formulas with the indices 
a,b,.... To simplify our formulas, we assume symmetry between the two 
neurons: Ru = R22, which in our case is very well satisfied - see Figure 6. 

The 4-point functions are now represented as 

'^1112 — [^^iuR^2 + ^^iaR\2 + ^ucxRl2]^lll2 + B^^jt^ 

'^1122 — [<Jt^u0al3 + ^12 ^12 + ^12 ^12 J^1122 + -Dll22 [^.60) 

'^Vni = [-^12^0/3 + Rl2^ul3 + Rl2 ^ua] ^1222 + B'^222' 
"^2222^ — l^fiu^a/S + ^fia^ul3 + <^/x/3<^i/o] ^2222 + -5222^- 

The intermediate steps 1 and 2 leading to equation 4.30 now increase, 
since we have to express several 4-point functions in terms of 2-point func- 
tions, not all of them being diagonal. In the particular case of the two HI 
neurons though, we may further simplify this system, neglecting R12. Its 
effect^ is very small indeed, since for rotational stimuli the action of the two 
neurons is complementary: an exciting stimulus for one neuron is inhibiting 
for the other - see Figure 6. Although for translational stimuli both neurons 
fire nearly synchronously, the dominant peak near r = in R12 is absent, 
since synchrony is not exact. In the following we therefore neglect K12. As 
can be seen in Figure 7, K12 is only ~ K22/5. Since the contributions of Ku 
and ^^22 are already small, Ki2^s 1 % effect can be safely neglected for both 
types of stimuli. 

Our equations now decouple and we get two sets identical to equations 
4.27, one for each neuron. 



^Thc effect of R^^ may be included perturbatively- see section 6. 
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5 Reconstructing the fly's stimulus and mea- 
suring its quality 

To test the quality of our reconstructions, we use the data with rj ~ 0.8, r = 
10 milhseconds and (r) ~ 80 spikes sec~^. 




32 48 64 16 32 48 

(bins) 



Figure 7: Second order kernel i^22(^i,^2) and upscaled version of -^'12(^1, ^2) for 
= 64. (A) 5 * Ki2 , (B) A'22. Notice i^22//'^i2 ~ 5. 

We select a representative sample, one second long, of the experiment, in 
order to give a visual display of the reconstruction. In Figure 8 we show the 
first order reconstruction of the original stimulus using Kl and K2 and the 
second order reconstruction, where the effect of Ku and i^22 is added - with 
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and without the Gl-approximation. We conclude: 

• Reconstructions using the experimental 4-point functions are very sim- 
ilar to their Gl-approximation. 

• The reconstruction procedure is unable to reproduce the fast stimulus 
variations at the 2 milliseconds time scale. It is also clear that still 
higher order terms are not going to improve this deficiency. But the 
second order kernels always represent an improvement, since the black 
line in Figure 8 is always a better approximation to the stimulus than 
the blue one. 

• We observe a stimulus-to-spike delay time of t^ot ~ 20 bins. 



S(t) 




1750 1800 1850 1900 1950 2000 2050 2100 2150 2200 2250 

time (bins) 



Figure 8: Reconstructing the rotational stimulus with kernels Ki, K2 and Kn, 
K22, using the experimental 4-point function and the Gl-approximation. Black 
thin dashed line: S{t), input stimulus to be reconstructed, blue line: Si{t), 
reconstruction using only Ki and K2, black continuous line: Si{t) + S2{t), 
experimental second order reconstruction, gray dashed line: Si{t) + S2{t) : Gl, 
Gl-second order reconstruction. The x and • signs stand for the right and left 
spikes respectively. Observe a delay-time of about 20 bins. 



15 



Although visual appraisement of the reconstruction quality is an indis- 
pensable guide to our intuition, numerical measures are less subjective. We 
naturally use the x^^"* = (/ dt[se{t) — of Equation 2.4, since its min- 

imization was used to determine the kernels ki,Kj. The reconstruction im- 
provement due to second order kernel is reflected in 

(2) _ (2) 

SX^'^ - (5-36) 
Xi 

where xf^ takes only first order terms into account - x'l^ = (/ dt[si(t) — 
s{t)Y), whereas second order terms are included in Xu — il dt[si{t) +S2{t) — 
s{t)f). 6x^^^ is posit ive, but small of ~ 8%. The chi-squared difference 
between the experimental and Gl-reconstructions is only of ~ 0.5 %. 

Although the x^^^-improvement is small, second order terms are a impor- 
tant at specific stimulus-dependent instants. In order to assess the relevance 
of these, we measure local chi-squares, defined as: 

rt+AT 

X?(t,AT)= / dt{{s,{t) - (5.37) 

Jt-AT 

and 

rt+AT 

Xl2it)= / dt{{si{t) + S2{t) - s{t))\) (5.38) 

Jt-AT 

for t = T2, where T2 are instants when Xui^) least as important as Xi(^)- 
If N2 is the number of such windows of size AT and Nt the duration of the 
experiment in bins divided by the window-size in bins, we plot in Figure 
9 the fraction of the stimulus-dependent instants vs. Xi/Xi2- Although 
this fraction vanishes as we require the importance of second order terms to 
increase, they still make a sizable contribution. Unfortunately just looking 
at the mean stimulus around T2 does not provide any insight and a more 
detailed analysis will be needed to reveal features, which might be relevant 
at these particular instants. 

Here we only follow (Rieke et al., 1997) and separate systematic from 
random errors, decomposing the estimate Se('^) into a frequency- dependent 
gain g{uj) and an effective noise neff{uj) referred to the input: 

Se{u;) = g{Ld)[s{uj) + ^^//(u;)]. (5.39) 
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Figure 9: N2/NT versus X1/X12 fo'' experimental and Gl-reconstruction. We 
find the instants where X1/X12 assumes a particular value > 1, when computed 
in windows of size AT = 64 bins. N2 is the number of these windows, whereas 
Nt is the duration of the experiment in bins divided by the window size. 



Around T2, we observe an overall improvement of 20% in g{uj). A further 
indication, that second order contributions, although drowned in averages 
over the whole experiment, may nevertheless have crucial importance in im- 
proving the code at specific moments. 

Finally we discuss the reconstruction of translational stimuli. Although 
in real life there is a continuous intermingling of rotational and translational 
motion, for a start we have considered this artificial separation of stimuli. 
Thus we have computed all averages (■) also for the translational setup. The 
kernels Ka-, Kab are similar to the rotational ones, but there is a sign change. 
Whereas for rotational stimuli Ki ~ —K21 Kn ~ —K22, for the translational 
case we have 

The reconstructions shown in Figure 10 are worse than the rotational ones. 
For positive stimuli, corresponding to unrealistic backward motion of the 
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Figure 10: Reconstructing tlie translational stimulus witli kernels Ki, K2 and 
Kii, K22, using the experimental 4-point function and the Gl-approximation. 
Black thin dashed line: S{t), input stimulus to be reconstructed, blue line: Si{t), 
reconstruction using only Ki and K2, black continuous line: Si{t) + S2{t), 
experimental second order reconstruction, gray dashed line: Si{t) + S2{t) : Gl, 
Gl-second order reconstruction. The x and • signs stand for the right and left 
spikes respectively. Observe a delay-time of about 25 bins. 

fly, both neurons fire vigorously, whereas in the opposite case none does. 
Interestingly, the delay-time is now tfrans ~ 25 bins, about 5 bins larger 
than trot'- inspite of their mutual inhibition, the neurons manage to fire, 
albeit a little bit retarded. The Gl-representation works equally well for this 
case. It would be interesting to subject the fly to a more realistic mixture 
of rotational and translational motion without separating the two and then 
compute correlation functions etc. We intend to come back to this issue in 
the future. 
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6 Gl-approximation in population coding: tam- 
ing the matrix explosion 



Although the spike generation process of the HI neurons is not Gaussian, 
the parametrization 3.20 is unexpectedly good. Actually we don't know 
how to judge from the spike interval distribution, whether this surprise will 
happen or not. In fact, the interval distribution of the spike times looks more 
nearly Poisson, instead of Gaussian. We remark, that independent increment 
probability distributions, whether they are Poisson or not, never do justice 
to correlated spike trains. On the other hand, if the 2-point function R{t) 
is to be a suitable building block to represent the 4-point function, then the 
parametrization, equation 3.20, is uniquely selected to be the most general 
one respecting the symmetry of -R^^^(l, 2, 3, 4). 

Since first order computations treat each neuron independently and do not 
take their mutual correlations into account, in the future one certainly would 
want to perform second order reconstructions to study the fly's visual system 
for more than two neurons. Our Gl-approximation makes these computations 
much more feasible. It should also work for correlation functions involving 
neurons not belonging to the fly's lobula plate. 

In order to apply our Gl-approximation, we imposed the requirement 
= i?22 and we neglected Ru- This limitation may be relaxed in the 
following way^. One could set R12 = and use a different set of functions for 
each neuron, diagonalizing thus all 2-point functions TZaa and compute the 
coefficients T>ab- Then reexpand all variables in terms of one set of functions 
only and apply the procedure, which led to equation 4.31 for Ru 7^ 0. If this 
does not lead to a closed set of equations, small effects may always be taken 
into account by a perturbative scheme to arbitrary order. In fact, suppose 
we have solved equation 4.27 for some representation of TV^^^ - e.g. as we 
did in section 3. Incorporating R12 7^ and/or i?22 7^ ^ij.^ will change the 
7^-matrix into: 

n' = n + 6n, (6.41) 

with 6TZ supposedly small. The new equations to be solved are: 

'5.T=E^''r/^^':^, (6.42) 

cd,al3 

^Hcre we only provided an outline, leaving a detailed analysis for a future publication. 
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where V = V + 6V and V satisfies the unprimed equations 4.27. Expand- 
ing both sides of equation 6.42 to first order in the corrections, we get the 
equations 

- E ^^r/ = E 5^^ct, (6.43) 

cd,af3 cd,al3 

to be solved for the unknowns 6V. {—6TZ ■ V) replaces the left-hand-side of 
equation 4.27 and couples the neurons. The right-hand-sides of the above 
equation and equation 4.27 have the same form and can therefore be solved 
in the same manner. 

The Gl- approximation could also be useful for other systems and this 
would be a considerable step forward in implementing coding involving a large 
population of neurons. One of the problems in second order reconstructions 
involving many neurons is the size-explosion of the 4-point function matrices 
alluded to at equation 2.18. If, e.g. we record from four neurons using 
128 bin-sized windows, the length of the matrices to be inverted would be 
~ 128^ X 2^ ~ 10^^. With our approximation the size of the linear system to 
be solved grows only linearly with the number of neurons. 

In order to use our approximation, one would have to check the win- 
dowsize independence of the parameters Aat... and Bab... for some subset of 
the complete matrix-indices, to convince oneself of the adequacy of the ap- 
proximation. Since in our case the matrices were still manageable, we could 
compute the experimental 4-point functions to verify this point, but this will 
in general not be possible. 

7 Materials and Methods 

Flies, immobilized with wax, viewed two Tektronix 608 Monitors Ml, M2, 
one for each eye, from a distance of 12cm, as depicted in Figure 1. The mon- 
itors were horizontally centered, such that the mean spiking rates of the two 
neurons, averaged over several minutes, were equal. They were positioned, 
such that a straight line connecting the most sensitive spot of the compound 
eye to the monitor was perpendicular to the monitor's screen. The light in- 
tensity corresponds roughly to that seen by a fly at dusk (Steveninck, Lewen, 
Strong, Koberle, & Bialek, 1997). The stimulus was a rigidly moving vertical 
bar pattern with horizontal velocity v(t). We discretise time in bins of 2 mil- 
liseconds, which is roughly the refractory period of the HI neurons. The fly 
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therefore saw a new frame on the monitor every 6t = 2 milhseconds, whose 
change in position 6x was given by Sx{t) = v(t)6t. 

The velocity v{t) was generated by an Ornstein-Uhlenbeck process with 
correlation times Tc = 0, 5 and 10 ms i.e. the stimulus was taken from a 
Gaussian distribution with correlation function C{t) = e~^l'^'^ . Experimental 
runs for each Tc lasted 45 minutes, consisting of 20 seconds long segments. In 
each segment, in the first 10 seconds the same stimulus was shown, whereas 
in the next 10 seconds the fly saw different stimuli. 

8 Summary 

The ability to reconstruct stimuli from the output of sensory neurons is a 
basic step in understanding how sensory systems operate. If intra- and inter- 
neuron correlations between the spikes emitted by neurons are to be taken 
into account, going beyond first order reconstructions is mandatory. In this 
case one has to face the size-explosion of higher order spike-spike correlation 
functions, the simplest being the 4-point correlation function necessary for a 
second order reconstruction. Our Gl-representation of the 4-point function 
in terms of 2-point functions tames this problem. If this representation holds 
more generally, the coding in large populations would become more feasible. 

For our case of the two HI neurons of the fly, correlations between them 
may be neglected, since they are only of ~ 1 %. We perform reconstructions 
using both the experimental and the Gl-approximation for the 4-point func- 
tions involved. Both are very similar, their chi-squared differing by 0.5 %. 
To implement the Gl-program for the two neurons, we found it convenient 
to expand our variables in terms of eigenf unctions of 2-point matrices. We 
propose a perturbative scheme in order to take the neglected correlations 
into account. 

We find that second order terms always improve the reconstruction, al- 
though measured by a chi-squared averaged over the whole experiment this 
improvement is only at the 8% level. Yet these terms can represent a 100 % 
improvement at special instants as measured by an instant dependent chi- 
squared. 



^"Although we show results only for = 10 ms our conclusions arc also valid for Tc = 0, 5 
ms. 
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